! gfortran -m64 -msse -msse2 -msse3 -m3dnow -mfpmath=sse \
!	-fdefault-real-8 -fdefault-integer-8 -fopenmp \
!	-pipe -Wall fpt.f90 -o fpt
program fpt
  implicit none
  integer,parameter::narg=5
  integer,parameter::max_scale=16384
  integer,parameter::MAX_PATHLEN=1024
  character(MAX_PATHLEN) fin,fout,str
  integer n
  integer nx,ny,nxy
  integer mx,my,mz,mz0
  integer ix1,iy1,ix2,iy2
  integer iy2f,ix2f
  integer jx,jy,jz
  real D,z,f,mv
  integer*2,allocatable::hist(:,:,:,:)
  real,allocatable::img(:,:,:)
  integer*2,allocatable::img2(:,:,:)
  integer istat

  nx=0; ny=0; D=0
  n=command_argument_count()
  if(n<narg) stop "gimme fin fout nx ny D"
  call get_command_argument(1,fin)
  call get_command_argument(2,fout)
  call get_command_argument(3,str)
  nx=a2i(str)
  if(nx<=0) stop "*** Invalid nx"
  call get_command_argument(4,str)
  ny=a2i(str)
  if(ny<=0) stop "*** Invalid ny"
  call get_command_argument(5,str)
  D=a2i(str)
  if(D<=0) stop "*** invalid D"

  nxy=nx*ny
  allocate(hist(nx,ny,nx,ny),stat=istat)
  if(istat/=0) stop "*** Error during allocate"
  
  open(unit=1,file=fin,iostat=istat,status="old",access="stream")
  if(istat/=0) stop "*** Error during open"

  read(1,iostat=istat) hist
  if(istat/=0) stop "*** Error during read"
  close(1)
  
  mx=nx*2; my=ny*2; mz=D
  if(iand(mz,1)==0) mz=mz-1 ! mz always odd
  mz0=mz/2+1
 
  allocate(img(mx,my,mz))
  img=0
  
  !$omp parallel do schedule(dynamic) private(jz)
  do jz=1,mz
     z=jz-mz0
     f=1.0/2.0-z/D
     do iy2=1,ny
        iy2f=iy2*(1.0-f)
        do ix2=1,nx
           ix2f=ix2*(1.0-f)
           do iy1=1,ny
              jy=iy1*f+iy2f
              do ix1=1,nx
                 jx=ix1*f+ix2f
                 img(jx,jy,jz)=img(jx,jy,jz)+1
              end do
           end do
        end do
     end do
  end do
  !$omp end parallel do

  allocate(img2(mx,my,mz))
  mv=maxval(img)
  if(mv>max_scale) img=img/mv*max_scale
  img2=img

  open(unit=1,file=fout,status="replace",access="stream",iostat=istat)
  if(istat/=0) stop "*** Error during open"
  write(1,iostat=istat) img2
  if(istat/=0) stop "*** Error during write"
  close(1)

  contains
    
    integer function a2i(s,ok)
      character*(*),intent(in)::s
      logical,intent(out),optional::ok
      integer istat
      integer x
      read(s,*,iostat=istat) x
      a2i=x
      if(istat==0) then
         if(present(ok)) ok=.true.
      else
         if(present(ok)) ok=.false.
      end if
    end function a2i

end program fpt
